\defmodule{Rank1Lattice}

This class implements point sets specified by integration
lattices of rank 1. They are defined as follows \cite{vSLO94a}.
One selects an arbitrary positive integer $n$ and a $s$-dimensional
integer vector $(a_0,\dots,a_{s-1})$.
[Usually, $a_0=1$ and $0 \le a_j < n$ for each $j$;
when the $a_j$ are outside the interval $[0,n)$, then we replace  $a_j$ by
($a_j \bmod n$) in all calculations.] The points are defined by
\eq
  \mathbf{u}_i = (i/n)(a_0, a_1, \ldots, a_{s-1}) \bmod 1
\endeq
for $i=0,\dots,n-1$.
These $n$ points are distinct provided that $n$ and the $a_j$'s have
no common factor.


\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        Rank1Lattice
 * Description:  Rank-1 lattice
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.hups;\begin{hide}
import umontreal.iro.lecuyer.util.PrintfFormat;
import umontreal.iro.lecuyer.rng.RandomStream;
\end{hide}

public class Rank1Lattice extends PointSet \begin{hide} {

   protected int[] genAs;          // Lattice generator:  a[i]
   protected double[] v;           // Lattice vector:  v[i] = a[i]/n
   protected double normFactor;    // 1/n.
   protected double[] shift;       // Random shift, initially null.

   private void initN (int n) {
      numPoints = n;
      normFactor = 1.0 / (double) n;
      for (int j = 0; j < dim; j++) {
         int amod = (genAs[j] % n) + (genAs[j] < 0 ? n : 0);
         v[j] = normFactor * amod;
      }
   }
\end{hide}
\end{code}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Constructor}
\begin{code}

   public Rank1Lattice (int n, int[] a, int s) \begin{hide} {
      dim = s;
      v = new double[s];
      genAs = new int[s];
      for (int j = 0; j < s; j++) {
         genAs[j] = a[j];
      }
      initN (n);
   }\end{hide}
\end{code}
 \begin{tabb}
   Instantiates a \class{Rank1Lattice}{} with $n$ points and lattice
   vector $a$ of dimension $s$.
 \end{tabb}
\begin{htmlonly}
   \param{n}{there are n points}
   \param{a}{the lattice vector}
   \param{s}{dimension of the lattice vector a}
\end{htmlonly}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Methods}
\begin{code}

   public void setNumPoints (int n) \begin{hide} {
      initN(n);
   }\end{hide}
\end{code}
\begin{tabb}
  Resets the number of points of the lattice to $n$. The dimension  $s$ and
  the $a_j$ are unchanged.
\end{tabb}
\begin{code}

   public int[] getAs() \begin{hide} {
      return genAs;
   }\end{hide}
\end{code}
\begin{tabb}
Returns the generator $a_j$ of the lattice. (The original ones before they are
reset to $a_j \bmod n$). Its components
  are returned as \texttt{a[$j$]}, for $j = 0, 1, \ldots, (s-1)$.
\end{tabb}
\begin{code}

   public void addRandomShift (int d1, int d2, RandomStream stream) \begin{hide} {
      if (null == stream)
         throw new IllegalArgumentException (
              PrintfFormat.NEWLINE +
                  "   Calling addRandomShift with null stream");
      if (0 == d2)
         d2 = Math.max (1, dim);
      if (shift == null) {
         shift = new double[d2];
         capacityShift = d2;
      } else if (d2 > capacityShift) {
         int d3 = Math.max (4, capacityShift);
         while (d2 > d3)
            d3 *= 2;
         double[] temp = new double[d3];
         capacityShift = d3;
         for (int i = 0; i < d1; i++)
            temp[i] = shift[i];
         shift = temp;
      }
      dimShift = d2;
      for (int i = d1; i < d2; i++)
         shift[i] = stream.nextDouble ();
      shiftStream = stream;
   }\end{hide}
\end{code}
\begin{tabb}  Adds a random shift to all the points of the point set,
  using stream \texttt{stream} to generate the random numbers.
  For each coordinate $j$ from \texttt{d1} to \texttt{d2-1},
  the shift $d_{j}$ is generated uniformly over $[0, 1)$ and added modulo $1$ to
  all the coordinates of all the points.
%  After adding a digital shift, all iterators must be reconstructed or
%  reset to zero.
\end{tabb}
\begin{htmlonly}
   \param{d1}{lower dimension of shift}
   \param{d2}{upper dimension of shift is d2 - 1}
   \param{stream}{random number stream used to generate uniforms}
\end{htmlonly}
\begin{code}

   public void clearRandomShift() \begin{hide} {
      super.clearRandomShift();
      shift = null;
   }\end{hide}
\end{code}
\begin{tabb}  Clears the random shift.
\end{tabb}
\begin{code}
 \begin{hide}

   public String toString() {
      StringBuffer sb = new StringBuffer ("Rank1Lattice:" +
                                           PrintfFormat.NEWLINE);
      sb.append (super.toString());
      return sb.toString();
   }


   public double getCoordinate (int i, int j) {
      double x = (v[j] * i) % 1.0;
      if (shift != null) {
         if (j >= dimShift)   // Extend the shift.
            addRandomShift (dimShift, j + 1, shiftStream);
         x += shift[j];
         if (x >= 1.0)
            x -= 1.0;
         if (x <= 0.0)
            x = EpsilonHalf;  // avoid x = 0
       }
      return x;
   }


   // Recursive method that computes a^e mod m.
   protected long modPower (long a, int e, int m) {
      // If parameters a and m == numPoints could be omitted, then
      // the routine would run much faster due to reduced stack usage.
      // Note that a can be larger than m, e.g. in lattice sequences !

      if (e == 0)
         return 1;
      else if (e == 1)
         return a % m;
      else if ((e & 1) == 1)
         return (a * modPower(a, e - 1, m)) % m;
      else {
         long p = modPower(a, e / 2, m);
         return (p * p) % m;
      }
   }

   protected double radicalInverse (int base, int i) {
      double digit, radical, inverse;
      digit = radical = 1.0 / (double) base;
      for (inverse = 0.0; i > 0; i /= base) {
         inverse += digit * (double) (i % base);
         digit *= radical;
      }
      return inverse;
   }

   public PointSetIterator iterator() {
      return new Rank1LatticeIterator();
   }

// ************************************************************************

   protected class Rank1LatticeIterator extends PointSet.DefaultPointSetIterator
   {
      public double nextCoordinate() {
         // I tried with long's and with double's. The double version is
         // 4.5 times faster than the long version.
         if (curPointIndex >= numPoints || curCoordIndex >= dim)
            outOfBounds();
//      return (curPointIndex * v[curCoordIndex++]) % 1.0;
         double x = (curPointIndex * v[curCoordIndex]) % 1.0;
         if (shift != null) {
             if (curCoordIndex >= dimShift)   // Extend the shift.
                addRandomShift (dimShift, curCoordIndex + 1, shiftStream);
             x += shift[curCoordIndex];
             if (x >= 1.0)
                x -= 1.0;
             if (x <= 0.0)
                x = EpsilonHalf;  // avoid x = 0
         }
         curCoordIndex++;
         return x;
      }
   }
}\end{hide}
\end{code}
